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FOREWORD 
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tation of an Internal Flow Field in a Scramjet Engine," for the period ended 
August 31, 1986. The work was supported by the NASA/Langley Research Center 
(Computational Methods Branch of the High-Speed Aerodynamics Division) 
through research grant NAG-1-423. The grant was monitored by Dr. Ajay Kumar 
and Mr. J. Philip Drummond of the High-Speed Aerodynamics Branch. 
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INVESTIGATION OF CHEMICALLY REACTING 
AND RADIATING SUPERSONIC INTERNAL FLOWS 

M. Mani* and S. N. Tiwari* 


ABSTRACT 

The two-dimensional spatially elliptic Navier-Stokes equations are used 
to investigate the chemically reacting and radiating supersonic flow of the 
hydrogen-air system between two parallel plates and in a channel with a ten- 
degree compression-expansion ramp at the lower boundary. The explicit unsplit 
finite-difference technique of MacCormack is used to advance the governing 
equations in time until convergence is achieved. The chemistry source term in 
the species equation is treated implicitly to alleviate the stiffness 
associated with fast reactions. The tangent slab approximation is employed in 
the radiative flux formulation. Both pseudo-gray and nongray models are used 
to represent the absorption-emission characteristics of the participating 
species. Results obtained for specific conditions indicate that the radiative 
interaction can have a significant Influence on the flow field. 


*Graduate Research Assistant, Department of Mechanical Engineering and 
Mechanics, Old Dominion University, Norfolk, Virginia 23508. 

^Eminent Professor, Department of Mechanical Engineering and Mechanics, 
Old Dominion University, Norfolk, Virginia 23508. 
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I 

1. INTRODUCTION 

In the last several years there has been a great deal of research toward 

i 

development of a hypersonic transatmospheric vehicle. At the NASA Langley 

| Research Center the hydrogen-fueled scramjet (supersonic combustion ramjet) 

j 

| engine has been a strong candidate for propelling such a vehicle. For a 

better understanding of the complex flow field in different regions of the 
engine, both experimental and computational techniques are being employed. 
Several computer programs Cl— 43* have been developed to gain more insight into 
the problem involving the flow in the various sections of the scramjet module. 

The purpose of this study is to investigate the effect of radiative heat 
transfer in chemically reacting supersonic flow in the scramjet combustor. 
The combustion of hydrogen and air results in absorbing-emitting gases such as 
water vapor and hydroxyl radical. The presence of such gases makes the study 
of radiative heat transfer in chemically reacting flows an important issue. 

There are several models available in the literature to represent the 
absorption-emission characteristics of molecular gases [5-11]. Some specific 
applications of the models to flow and combustion related problems are 
available in [8, 9, 12-14], Both gray and nongray gas models are used in this 
study. The gray gas model is less accurate but is much more efficient for 
parametric studies. For this model, the radiative heat flux is Independent of 
the wavelength and the governing equations can be expressed in terms of second 
order nonhomogenious ordinary differential equations (ODE). This system of 
ODE's forms a tridiagonal matrix which is solved by the Thomas algorithm. 

The flow field in the combustor is represented by the Navier-Stokes 


*The numbers in brackets indicate references. 
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equations and the appropriate number of species continuity equations [1-4]. 
Incorporation of the finite rate chemistry model into the fluid dynamic 
equations creates a set of stiff differential equations. The stiffness is due 
to a disparity in the time scales of the governing equations. In the time 
accurate solution after the fast transients have decayed and the solutions are 
changing slowly, taking a larger time step is necessary for efficiency 
purposes. But explicit methods still require small time steps to maintain 
stability. An eigenvalue problem associated with stiff ODE has been solved to 
express this point clearly in [2]. One way around this problem is to use a 
fully Implicit methods. This method, however, requires the inversion of a 
block multi-diagonal system of algebraic equations. This is difficult to 
implement to take full advantage of vector processing computers such as VPS- 
32. The use of a semi-implicit technique, suggested by several investigators 
[15-17], provides an alternative to the above problems. This technique treats 
the source term (which is the cause of the stiffness) implicitly, and solves 
the rest of equations explicitly. 

2. GENERAL FORMULATION 

A brief discussion is presented on various components of the scramjet 
engine.' Special attention is directed in discussion of the basic equations 
that are applicable in analyzing the flow field in different parts of the 
engine. The relations for the thermodynamic and chemistry models are also 
provided in this section. 

2.1 Physical System and Model 

As mentioned in the introduction, the scramjet engine has been a 
candidate for propelling the hypersonic vehicle. In Fig. 2.1, various air 
breathing and rocket propulsion alternatives are shown [18, 19]. For Mach 
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numbers zero to three turbojet air breathing systems have the highest 
performance. Above Mach number of three, turbine inlet temperatures constrain 
performance, and then the ramjet becomes more attractive. At about Mach 
number of six, the performance of the ramjet is greatly reduced. This is due 
to dissociation of the reaction product which is caused by slowing the 
supersonic flow to subsonic flow. Therefore, it is more efficient to allow 
the engine internal flow to remain at supersonic speed. Thus for Mach number 
of six and beyond, the fixed geometry scramjet is clearly superior for 
propelling a vehicle at hypersonic speed. Hydrogen has been selected as the 
fuel for the scramjet due to its capability of cooling the engine and the 
airframe and also because of its high impulse level. 

The scramjet engine is made up of several identical modules and it is 
installed underneath the aircraft as shown in Fig. 2.2. As part of the engine 
design concept, the forebody of the aircraft acts as an inlet (airframe 
integrated [19]) for precompression and the afterbody as a nozzle for post 
expansion. Since the vehicle compresses the airflow in the vertical 
direction, the module size walls are made wedge shaped to compress the flow 
horizontally. This tends to minimize the flow distortion. The diamond shape 
strut is located at the minimum cross sectional area to complete the 
compression and it also provides the fuel injection. Each module is made up 
of inlet, combustor and nozzle regions. 

The inlet region starts with the forebody of the vehicle and ends up with 
the minimum cross sectional area of each module. In the first part, the air 
flow is compressed by the oblique shock generated from the forebody before it 
enters the engine. For numerical solution, the flow is best represented by 
the Navier-Stokes equation in the actual inlet area of the engine. Using the 
Euler equation away from the wall and the boundary layer equation near the 


4 




Fig. 2.2 Airfram-integrated supersonic combustion ramjet 
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wall region can be complicated by the oblique shock interaction with the 
boundary layer. This can cause the flow separation which means flow can not 

be represented accurately by these equations. The chemistry is frozen in this 

% 

region (t ch >>T f j Three-dimensional Navier-Stokes equations have been 

employed by Kumar [3] to investigate the flow field in this region with a 
reasonable success. Chltsomboon et al. [4] have employed the parabolized 

Navier-Stokes equations with limited success. 

The combustor region is by far the most complex part of the scramjet 

engine. As a result, a great deal of research is directed toward better 

understanding of the combustor flow field. The flow in this region is usually 
supersonic but it can be both supersonic and subsonic (Fig. 2.3). The fluid 
dynamics becomes complicated by the fuel injection, flame holding, chemistry, 
radiation and turbulence. The flow field in this region is represented by the 
elliptic Navier-Stokes (including turbulence, chemistry and radiation) 
equations. In the farfield (downstream of the fuel injection strut where the 
flow is only supersonic) the flow can be represented by the parabolized 

Navier-Stokes equations [4], 

The nozzle and subsurface of the afterbody provide about fifty percent of 
the thrust at Mach number six [19]. The flow through the nozzle is supersonic 
and the chemistry is frozen. However the combustor exit flow consists of 
multi components of the reacting species, multiple shock, and 3-D viscous 
effects. The flow can be represented by the parabolized Navier-Stokes 
equations. 


2.2 Basic Governing Equations 

The two-dimensional Navier-Stokes and species continuity equations are 


represented in physical domain by 




Fig* 2.3 Flow field near an Injec 
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5U dF dG „ _ n 
FE + 5x + 3y + H ’ 0 


( 2 . 1 ) 


where vectors U, F, G and H are expressed as 


U = 


P 

pu 

pv 

pE 

Pf, 


F = 


G = 


pu 

PU 2 + T, 


XX 

puv + X xy 

(PE ♦ t xx ) u + x xy v + q cx 


5 f . 

pu, j ' ^ sr 


pv 

puv 2 + T yx 

pv +T yy 

(pE + x )v+t u+q +q„ 
yy ftf xy M cy M Ry 

pvf j " pD W 


H = 


0 

0 

0 

0 


The viscous stress tensors in the F and G terms of Eq. (2.1) are given as 


*xx = p -< + -57 ) ~ 2pdU 


dx 


(2.2a) 


_ , ,/6u . dv» 

T xy ” ' p( dy + 


(2.2b) 


t = p - x(2!U& - 2p|l 
yy TSx ^y K 75y 


(2.2c) 


The quantities q cx and q C y in the F and G terms are the components of the 
conduction heat flux and are expressed as 


<U " - k 


cx 


-rrr “ pD E C(bf^/9x)h.] 
OX j =1 J J 


dT 


m 


Icy * * k W ' p0 jfj Cd V 8y) V 


(2.3a) 


(2.3b) 


"j ■ h j ♦ I C P ( dT ; T o - 0 K 

* A J 


where 
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It should be noted here that D represents the effective binary diffusion 
coefficient and is used for all species. Assuming that the Lewis number is 
unity, Eqs. (3) reduce to (see Appendix A) 

vfl < 2 -*» 

(2 - 4bl 

where 

e = h - P/p 

The molecular viscosity \l is assumed to be temperature dependent and it is 
evaluated for individual species from the Sutherland's formula 

3/2 T + S 

p = _ T + s * 2,5 * 

o 

where p Q and T 0 are reference values and S is the Sutherland constant. The 
total internal energy E in Eq. (2.1) is given by 

2 2m 

E = - P/p + - — t v + 2 h. f i (2.6) 

* ' i*l 1 1 

Specific relations are needed for the chemistry and thermodynamic models 
and for the radiative transport. The chemistry and thermodynamic models are 
discussed briefly in the following sections. The formulations for the 
radiative transport are presented in Sec. 3. 

2.3 Chemistry Model 

For the numerical solution of reacting flows, a chemistry model is needed 
to represent the combustion process. The chemistry model used in this study 
is the two-step global finite rate hydrogen-air combustion model developed by 


9 


Rogers and Chlnitz [20]. This chemistry model was deduced from a 28 reaction 
model and it is adequate for temperatures between 1,000 K and 2,000 K and for 
equivalence ratios between 0.2 and 2.0. In the first step, hydrogen and air 
react and produce hydroxyl radical and in the second step, hydroxyl radical 
and hydrogen react to produce water vapor. The reactions are expressed as 


\ 

H 2 + 0 2 — i. 2 


OH 


'1 


20H + H 2 



(2.7) 


( 2 . 8 ) 


where K, 


and 


(. represent the forward and reverse reaction rate 
■i b i 

constants respectively. The relations for K* are obtained from an Arhenius 

T i 

equation as 


N 


i 


( f = A^iOT ' exp(-E i /RT) 


(2.9) 


The values of the parameters (4»), N^, and Ej in Eq. (2.14) are 

A 1 (4>) = (8.917 4> + 31.433/4. 

- 28.95 )x 10 47 ~ cm 3 /g-mol-S 
E x * 4865 Cal /mol; ^ * -10. 

A 2 (4>) = (2. + 1.333/4> 

- .8334>)xl0^ 4 ~ cm^/mol 2 - s 
E 2 = 42500 Cal /mol; N £ = -13 
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where $ is the fuel -air ratio. The reverse rate constants can be evaluated by 


! »i ' V 


( 2 . 10 ) 


where 


K x = 26.164 exp (-8992/T ) 


,-6 


K 2 = (2.682x10 °) (T) exp(69415/T) 

Knowing the reaction rates (K. and K k J, the production of species can be 

T i D i 

evaluated from the law of mass action. Consider the general chemical reaction 


N K f, N 

jfj v i.i c j ^ A v j!i c i , ' 1 - 2 m 

\ 


(2.11) 


where . and v* ‘ represent the stoichiometric coefficents of the reactants 
and products respectively. The law of mass action states that the rate of 
change of concentration of species j by reaction i is given by [21] 


, N Vjj N v'A 

(£..) - (vil - v’.. ) [K f n c . J1 - k. n c. j1 ] 
J jl Jl f i j= i J b i j=1 J 


( 2 . 12 ) 


The net rate of production of species j in all reactions is given by 


m 


C. = L (£.). 

J i=1 J 1 


(2.13) 


where m is the number of reactions. Finally, the chemistry source terms, on a 
mass basis, are found by multiplying the molar changes and corresponding 
molecular weight 


(2.14) 


11 


By applying the law of mass action to the global model, the chemistry source 
terms of the four species are obtained as 


^0 2 * " K fl C H 2 C 0 2 + K b x C 0H 


(2.15) 


^H 2 0 = 2(K f2 C 0H C H 2 ” K b 2 C H 2 0 } 


(2.16) 


^ H 2 = C 0 2 " 1/2 


(2.17) 


^OH = “ (2 ^0 2 + ^H 2 0 } 


(2.18) 


2.4 Thermodynamic Model 


The specific heat of individual species, C p , assumed to be a linear 
function of temperature, i.e.. 



* a. T + 

J 


(2.19) 


where aj and bj are constants which are obtained by curve fitting the 
thermochemical data of Ref. 22. The numerical values of these constants are 
given in Table 1. The specific heat of the mixture are computed by summing 
specific heat of individual species weighted by species mass fraction 

m 


C = Z C f, (2.20) 

p i=i p j J 


The static enthalpy of the mixture can be expressed as 

m T 

h = Z lh°. + / C dT] f. (2.21) 

j=l J T 0 P 1 3 

The total enthalpy can now be evaluated as 



12 


H * h + 0.5 (u 2 + v 2 ) (2.22) 

Combining Eqs. (2.21) and (2.22), the total enthalpy is expressed as 

2 

m a . T 4 o 

H = Z [h. + -V- + b. T] f . + 0.5 ( u + v*) (2.23) 

i=l 1 £ 11 

where h° is the sensible enthalpy of individual species at a reference 
temperature (T Q = 0 K). The gas constant for the mixture also is evaluated 
by a mass weighted summation over all the species as 

m 

R = Z f. R. (2.24) 

j=l J J 

The equation of state for the mixture of the gases, therefore can be written 
as 


P = PRT 


(2.25) 
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Table 2.1 Numerical Yalues of Various Constants 


Species 

H°(J/kg) 

a 

b 

o 2 

-271267.025 

0.119845 

947.937 

h 2 o 

-13972530.24 

0.43116 

1857.904 

h 2 

-4200188.095 

2.0596 

12867.46 

OH 

+1772591.157 

0.16564 

1672.813 

n 2 

-309483.98 

0.10354 

1048.389 
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3. RADIATION TRANSPORT MODELS 

In order to include the effects of radiative interaction in a physical 
problem, it is essential to accurately model the absorption-emission 
characteristics of participating species and provide a correct formulation of 
the radiative transfer processes. These are discussed briefly in this 
section. 


3.1 Radiation Absorption Models 

Many models are available in the literature to represent the absorption- 
emission characteristics of molecular species; a review of important models is 
available in [11]. Perhaps the simplest model is the gray gas model where the 
absorption coefficient is assumed to be independent of the wavelength. Many 
nongray models are also available in the literature. Both gray and nongray 
absorption models are discussed here briefly. 


3.1.1 Gray Gas Models 

In gray model, it is assumed that the absorption coefficient is 
independent of wavelength. This is rarely a physically realistic approxi- 
mation but it serves as an initial step for studying the effect of radiative 
heat transfer. The absorption coefficient for gray gas is evaluated by 
employing the Planck mean absorption coefficient as follows 


< 

P 


k e. (T) du 
a> bu> 
o 


e. (T) 

D 


(3.1) 


By assuming that within a band the Planck function does not vary significantly 
with the wave number and evaluating its value at the band center, the relation 
for <p for a single-band gas can be written as 
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where u 

o 

given by 


c »“„ <n 

<_ * * — / k dcj 

^ a T (y) Aw w 


(3.2) 


represents the band center. For a multiband gaseous system, < p is 


>L <T1 

* [- 4 —/ 


n 


T J K d<*»] 

k=l a r (y) Ao) w 


(3.3) 


where u> k represents the band center of the kth band of a particular species. 
For a specific band of a given gas, the Integrated band intensity S k is 
defined as 

s “TT ! da> (3 * 4) 

J Au 

Substituting Eq. (3.4) into Eq. (3.3), x p is expressed as 


where 


K P = 


P. n 

— i — 2 e b«* (T) s k (T) 
a r(y) k=l k 


r 3 

e bio |c = exp lc 2 c^/tj -i 


(3.5) 


The Pj is species partial pressure and Ci and C 2 are constants. Note that 
Kp is a function of temperature and species partial pressure. 

3.1.2 Nongray Gas Models 

Important nongray models available in the literature are as follow: 

1. Line Models: 

(a) Lorentz 

(b) Doppler 



16 


(c) Loren tz-Doppler (Voigt) 

2. Narrow Band Models 

(a) Elsasser 

(b) Statistical 

(c) Random-El sasser 

(d) Quasi-Random 

3. Wide Band Models 

(a) Box or Coffin 

(b) Modified box 

(c) Exponential 

(d) Axial 

The relative importance and range of applicability of these models are 
discussed in [11]. In the moderate temperature range (500-5000K) use of the 
wide band models and correlations provide sufficient accuracies. These models 
render significant computational efficiency over the line by line or narrow 
band models. 

The expression for the total band absorptance is given as 

00 _i 

A(y) = / [l-exp(-< y)] dw ~ cm (3.6a) 

o 

where both, w and ic^ have units of cm"*. Differentiation of Eq. (3.6a) gives 

00 _ p 

A'(y) = / x exp (-< v) du ~ cm (3.6b) 

o 

A"(y) * J - k* exp (-x^y) du> ~ cm* 3 (3.6c) 

The radiative flux term usually involves multiple integrals even for the 
simple geometries. As result numerical calculation of radiative flux for 
energy transfer between two parallel plates becomes very time consuming. 
Therefore it is desirable to replace the relation for the total band 
absorptance as given by Eq. (3.6a) with a continuous correlation [6-9]. 
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Numereous correlations are available in the literature for the wide band 
absorptance. The first correlation to satisfy the linear, square-root, and 
logrithmic limits of the wide band absorptance was proposed by Edward and 
Menard [23]. The most widely used correlation is the Tien and Lowder 
continuous correlation and this is given by [6] 

A = tn'Cuf(p) + 1.] (3.7) 

where 

f(p) * 2.94 [l-exp(-2.6p)] 

0 - B 2 P e 

The form of f(p) was chosen to give agreement with the correlation of Edward 
and Menard. This correlation is employed in this study for nongray gas 
formulation. 

3.2 Radiative Flux Equations 

The equations of radiative transport are expressed generally in integro- 
differential form; the integration involves both frequency spectrum and 
phyical coordinates. In many realistic three-dimensional physical problem, 
the complexity of the radiative transport equations can be overcome by 
introduction of the "trangent slab approximation." This approximation treats 
the gas layers as a one-dimensional slab in evaluation of the radiative 
flux. Tangent slab approximation is employed in this study. Therefore, the 
radiative transport is considered only in the normal direction of the flow. 
It should be pointed out that this approximation is not used for any other 
flow variables. 
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3.2.1 Basic Formulation 

The radative transport equations in the present study are obtained for a 
gas confined between two parallel plates (Fig. 5.1). For one-dimensional 
absorbing-emitting medium with diffuse boundaries the general equation for the 
radiative flux under the condition of thermodynamic equilibrium (LTE) is given 
by [5, 8, 24]. 

q RX = 2 B 1X E 3 (t X ) " 2 B 2X E 3 U 0X “ T X* 


X n / .\ r- / i \ j i a t OX 


+ 2n / * B x (t) E 2 (x x -t) dt - 2it J 0A> B x ( t) EoU-t ) dt 
0 T x 


(3.8) 


where 


V s \\ = 


The quantities B^ and B 2 ^ in Eq. (3.8) represent the spectral surface 
radios! ties, and for non-reflecting surfaces B lX =e lX =e lX e blX- For 
non-reflecting surfaces under the conditions of LTE the expression for the 
spectral radiative flux is expressed in terms of the wave number as (see 
Appendix B) 

= e lw ~ e 2u , + 2{ J o ^ e u> " e lJ E 2 (t u " t} dt 


'lu 2b) 
x 


- / °“ Ee (t) - e,J E, (t - x ) dt) 


(I) 


2u> 2 


b> 


(3.9) 


b) 


where E n ( t) is an exponential integral function defined by 


E n (t) = f e ' t/ ^ dp ; / E n (t) dt = - E n+1 (t) 


By employing the exponential kernal approximation [5] 
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E 2 (t) 


3 1 * E i (t) = i e 


_ 3 t 
9 7 1 


Eq. (3.9) is expressed as 






-7<V« 


dt 


where 


- / T °“ F 2u (t) e~ 7 (t ' T “’ dt) 

X 

0) 


F l„ (t) * « u <*> - *1„ i F 2u (t) - ejt) - e 


2d) 


Equation (3.11) is expressed in physical coordinate as 


V (y) = e l(, ' e 2 u + I f F lJ z) K * ex P C ' 1 K u <y - 2)]dz 

0 

• | / F 2w ( z) ex P [ ' 1 <J z -y)3dz 

By differentiating Eq. (3.12) with respect to y, one obtains 

- a q R Jy) q y 2 3 

— ar 1 — 3 7 f Q F iJ z) *1 ex Pt- 1 (y- z >]<* z 

//*.«*> 4 exp[-|K u (z-y)]dz 

- ! K <o [F U) (yl * F 2u < y >J 
The total radiative flux is given by 

oo 

V y) * I W y) d “ 

0 


(3.10) 


(3.11) 


(3.12) 


(3.13) 


(3.14a) 
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such that 

dq R (y) - dq Ru) d - 

3y = ! ~W d “ “ W t q Ro) du (3 ‘ 14b) 

0 0 

It should be noted here that for nongray gases, the divergence of 

radiative flux is used as a source term in the energy equation to avoid a 

scheme dependency in the computation. 

3.2.2 Gray Formulation 

In the previous section, it was observed that the radiative flux terms 
are represented by an integro-differential equation. Solving these equations 
are very time consuming, even with the vector processor such as VPS-32. 

Therefore, a pesudo-gray model is selected for efficient parametric studies. 
To express the radiative flux for a gray medium, one may assume that is 

now independent of the wave number. Therefore, Eq. (3.9) for a gray medium is 
written as 

T 

q R (x) = e 1 - e 2 + 2 / [e(t) - ej] E 2 (t - t)dt 

- 2 / 0 [e(t) - e 2 ] E 2 (t - t)dt (3.15) 

Upon substituting the exponential kernal approximation, Eq. (3.10), into Eq. 
(3.15), q R (-c ) is expressed as 

o x -4(x-t) 

q R (x) = e i " e 2 + 7 {/ Ce(t) - e^D e dt 


x -l(t-x) 

- / 0 Ce(t) - e 2 ] e 6 dt> (3.16) 
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To solve the radiative flux term for the gray medium, one can transform 
Eq. (3.15) into the physical coordinates and then apply any of the standard 
integration techniques to evaluate the radiative flux term. If this is 
desirable, then there is no need to make exponential kernal approximation and 
obtain Eq. (3.16). The main reason for employing the exponential kernal 
approximation in the gray formulation is to transform the integral equation 
into a differential equation. Solving a radiative flux equation in 
differential form is not only convenient but computationally efficient. For 
the present case, differentiating Eq. (3.16) twice by using the Leibnitz rule 
results in 


d q R (t) 

7 ? 


<|/ T Celt) 


-i (x-t) , i 

e^] e dt - 2 / 


< (t-T) 

[e(t) - e 2 ] e 6 dt> 


+ (e x - e 2 ) + 3 (3.17) 

A substitution of Eq. (3.17) into Eq. (3.16) gives a second order 
nonhomogenious ordinary differential equation as 


d q R (x) 


9 

? 


q R U) - 


de(x) 

dx 


(3.18) 


It should be pointed out that if Eq. (3.15) is differentiated and then Kernal 
approximation is applied the coefficient on the right side of Eq. (3.18) will 
be four. Employing the method of differential approximation instead of 
exponential kernal approximation, Eq. (3.15) may be recast as a differential 
equation which is of the same form as Eq. (3.18), except the coefficient of 
the second differential is 3/4 instead of unity [5], Equation (3.18) requires 
two boundary conditions, which for nonblack diffuse surfaces are given as (see 
Appendix C). 
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■ 7 1 q R 



1 dq R, 

J dx | x=0 



- e b (0) 


(3.19a) 


( ^"7 ) q R (T) |T=x 0 + ?7 t|t=t 0 * e b ( V ■ e b 2 (3 * 19b) 

Note that the right sides of Eqs. (3.19) differ from zero only if molecular 
conduction is neglected. In this study, other modes of energy transfer are 
included. Therefore, the terms on right sides reduce to zero. 


To evaluate the radiative flux, Eqs. (3.18) and (3.19) are transformed 
into the physical coordinates as 


d q R (y) q 

~ 7 d R (y) = ^ 

d(* p y) 2 4 R 


de (y) 
dU p y) 



(3.20a) 


(3.20b) 


(3.20c) 


The optical coordinate and thickness used in the above transformation are 
defined as 


* * * p y •• \ - « P L 

3.2.3 Nongray Formulation 

The radiative flux equation for nongray gases is obtained by substituting 
Eq. (3.13) into Eq. (3.14) as 

* I ! F i u (zl K l expC ' I \, <y ' z)]dz 
0 0 

+ f F 2(0 {z) K i expC " 7 K J Z ~ y )] d>do) 



23 


- 1 f *» tF l» ( »> + F 2, ( > )1(l * (3 - 21) 

0 

For a multi-band gaseous system, the above is expressed as 

W - " t A I <; * 'W 1 ’ 4 expC ' i V (y ' 2)1<lz 

J 1-1 Aa)^ Oil 1 

♦ ^ f 2», (z) \ expC ' i V | (z ' y)1 ' lz><i “i 

- 1 j, / t F to < y » + F 2«, < y >3 d “l < 3 - 22 > 

1 *1 Aw,j 1 1 

It should be pointed out that the following relation has been employed to 
obtain Eq. (3.22) 

I </ y F iJ 2 > expt'l Kjy - Zlldzjda., 

0 0 

N y ? 

= Z J {/ F, (z) tc exp[- i k (y - z)]dz)du. 

1=1 to, o J “i “l 2 “l 1 


where N represent the number of bands In a multiband system. Now, utilizing 
the definition of band absorptance and its derivatives as given in Eqs. (3.6), 
and evaluating the value of Planck function at the center of each band Eq. 
(3.22) is expressed as 


<Jq R (y) _ 3 n 

' 7 1=1 *"o1 


ay— = 1 ,?, <CF to_, (y) * V, (yn V d “i» 


01 


Aw.. 1 


N y 


+ T 2 if F, (z) a; 1 [| (y - z)] dz 
4 i=I o iu> oi 


+ f F 2oj ^ A i' ^7 ^ 2 “ y ) 3<i2> 

y oi 


(3.23) 
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Equation (3.23) is in proper form for obtaining the nongray solutions of 
molecular species. However, in order to be able to use the band model 
correlations, these equations must be transformed in terms of the correlation 
quantities. The correlation quantities and detail of transformations are 
given in Ref. 24. After the transformation, Eq. (3.23) is written as 

dq R (u) q N u. , 

-ar- ■ t £ A oi </ o F i Ul <“'> K c? »«i - «s 

♦ f 2 » ( (u '> K 4 <“i - 
U 1 1 

+ 1 ^ A o1%<“> +F 2 Ui <“” <3-24) 


Note that A“(u) express the second derivative of A(u) with respect to u and 

dq R _ dq R du _ r pc/ T w A -i dq R 

W "3IT w CPS(T)/A 0 ] IT 
By defining = £ and | Eq. (3.29) is expressed as 


dq R (y) _ 9 ! A o u o t ? c 
— 3y tt £ — * — U r 


T u 7 (z) ' ^7 T“ (y “ z) ]dz 

4 i=l L 2 o 1<0 i 1 2 L 


+ 1 *i' 4 <“i - “t> ]du i> 


+ ♦ F a s ( *» < 3 - 25 > 


It is often desirable and convenient to express the above equation in 
terms of A' rather than A 1 '. This is accomplished by integrating Eq. 
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(3.25) by parts. The details of the Integration by parts is given in Ref. 24 
and the result is given by 


dq R 3 " A °1 “°1 

W ' ? ,,j — r~ 



de (z) 
"i 
<3z 


A* 


[• 


3 % 
7~T 


(y - z)]dz 


- / A* [| -T-l (z - y)]dz} (3.26) 

Equation (3.26) and the Tien and Lowder correlation given by Eq. (3.7) can be 
used to evaluate the radiative flux. 


4. METHOD OF SOLUTION 

The grid generation and solution procedures for the governing equations 
are briefly discussed in this section. 

4.1 Grid Generation 

The grids are generated using an algebraic grid generation technique 
similar to the one used by Smith and Weigel [25]. From the computational 
point of view, it is desirable to have uniform rectangular grid enclosed in 
parallel piped, and the exterior of the parallel piped represent the physical 
boundaries. To have such grids, the body-fitted coordinate is transformed 
linearly from physical domain (x,y) to computational domain (5,n) as follow 


X x =7 (5,0) / 

y 2 * Y (5,0) ) 


Lower 

Boundary 


X 2 =7 (5,1) 


(4.1a) 


y 2 = Y (5,1) ) 


Upper 

Boundary 


(4.1b) 



Between the 
Boundaries 


(4.1e) 
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X = 7 (5,1) h + 7 (5,0) (1-n) I 
y * Y (5,1) n + y (5,0) (1-r,) J 


where 


0<5<1 ; 0 < n < 1 


The grid should be concentrated in the regions of high gradient to accurately 
capture the solution. Therefore more grids are required near the solid 
boundaries. The concentration of the grid in ti direction can be 
accomplished by 

_ (P y +1) - (p y -l) Exp[-C(rrl+a)/(l-a)] 

= I2a+1) U+ExpllClii-l+a)/ir-a)J} " (4 * 2) 


where 


P» + l 



If a is equal to zero (a=0) the compression takes place only near the 

lower wall (rpO), and if, a is set equal to one half (a*l/2), the 
compression takes place near both walls. The p y has a value between one and 
two, and as it gets closer to one, the grid becomes more concentrated near the 
walls. Employing this concentration, Eq. (4.1e) interms of ti is written as 


X = T (5,1) n + 1 U,0) (1-n) 

y = Y (5,1) n + Y (5.0) (1-q) (4.3) . 


where 


0 < T) < 1 


The grid used in this study is shown in Fig. 4.1. It should be noticed that 
the grid is concentrated in the normal direction and kept uniform in the flow 
direction. 
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Fig. 4.1 Grid Structure. 
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4.2 Solution of the Govenlng Equations 

The governing equations, Eqs. (2.1), are expressed in the computational 
domain as 


where 


5U + 3F 5G - _ Q 

'Sl + ^ + ^r\ +H 0 


U = UJ 

f - % - G \ 
A 

G = Gx^ - Fy^ 


H = HJ 


(4.4) 


0 - x 5 - \ \ 


Equation (4.4) is discretized temporally and written as 


u n+1 . u" - 4t [f + || + H n+1 ] 


i) n+1 = H n ♦ it 


■51 


H n+1 . H" + At » 

au Ir ^ 


(4.5) 

(4.6) 

(4.7) 


A substitution of Eq. (4.7) into Eq. (4.5) gives the temporally discrete 
equation in delta form as 


[I + At 51] AU n+1 = -At [|g + + H] n 


(4.8) 


where U n+1 - U n is expressed as AU n+1 , 51 is the Jacobian matrix of H and I is 


au 


the identity matrix. 


Once the temporal discretization used to construct Eq. (4.8) has been 
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performed, the resulting system is spatially differenced using the unsplit 
MacCormack predictor-corrector scheme [26]. This results in a spatially and 
temporally discrete, simultaneous system of equations at each grid point [16- 
17]. Each simultaneous system is solved using the Householder technique [27, 
28] in combination with the MacCormack technique which is then used to advance 
the equations in time. The modified MacCormack scheme then becomes 


[I + At ] 

au 1J 


dH-,n i *n+T 


AU 


ij 




(4.9a) 


r,n+T _ r.n . . r.n+1 
U u - U u + 4 u ij 


(4.9b) 


[ ! + At (|y AU- = -At[^ + ^ + H] ? : 


;n+l _ 


■5F . 5G . Cin+ 


(4.10a) 


U n+i 

U 1j 


r.n . n c r a • 

ij + LAU^j 




(4.10b) 


Equations (4.9) and (4.10) are used to advance the solution from time n 
to n+1 and this process is continued until a desired integration time has been 
reached. 


For time accurate solution, the computational time step. At, must satisfy 
the smallest time scales of the fluid and chemistry, i.e.. At = min 
(At f , At cfl ). If the steady state solution is sought (as in this study), it 
is possible to speed up the convergence by using a larger time scale 
(At^ * At^p|^) without introducing any instabilities in the solution. This is 
due to the so-called precondition matrix (left-hand side bracket in Eq. 
(4.8)), the purpose of which is to normalize the various time scales so that 
they are of the same order [17], 
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The radiative flux term Is evaluated for both gray and nongray gaseous 
systems. In the nongray gas formulation, the divergence of the radiative flux 
is evaluated using a central differencing scheme and is treated as radiative 
source term in the energy equation. Since the radiative flux term is in 
integro-differential form, unlike the other flux terms which are only in a 
differential form, it is uncoupled and treated separately. In the gray gas 
formulation, Eqs. (3.20) are discretized by central differencing, forming a 
tridiagonal matrix. This tridiagonal matrix can be solved efficiently by the 
Thomas algorithm. 

5. RESULTS AND DISCUSSION 

Based on the theory and computational procedure described in the previous 
sections, a computer code was developed to solve the two-dimensional Navier- 
Stokes equations for reacting and radiating supersonic laminar flows. Two 
different geometries are employed for various parametric studies. One is a 
channel with two parallel plates a distance L apart (Fig. 5.1); the other is a 
channel with a compression-expansion ramp at the lower boundary (Fig. 5.2). 
The freestream conditions at the inlet are obtained from Refs. 1-4. For the 
temperature range expected in the scramjet combustor, the radiating species 
that are important are OH and H 2 0. The spectral Information and correlation 
quantities needed for these species are obtained from Refs. 6-9 and 29. Both 
reacting and nonreacting flows are considered. Results for the nonreacting 
flows are presented in Figs. 5.3-5.10 and for the reacting flows results are 
presented in Figs. 5.11-5.16. 

For the parallel plate case (3 cm x 10 cm), the inflow conditions are 

P m = 1 atm, T^ = 1,700 K, = 3.0, and 0 = °’ 5, f 0 = 0,1 and 

f M = 0.4. Results for the radiative flux, as a function of the nondimen- 
n 2 
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sional location along the flow, are illustrated in Fig. 5.3 for various 

distances from the lower plate. It is noted that the radiation flux is 

approximately zero in the center of the plate (y = 1.5 cm) and is considerably 
higher towards the top and bottom plates. This, however, would be expected 
because of the symmetry of the problem and relatively higher temperature near 
the boundaries. The oscillations in results are due to the shock reflection 
from the boundaries. 

The results for radiative flux are illustrated in Figs. 5.4 and 5.5 as a 
function of the nondimensional y-coordinate. For P = 1 atm, the results 
presented in Fig. 5.4 for different water vapor concentrations indicate that 
the radiative interaction increases slowly with an increase in the amount of 
the gas. The results for 50* H 2 O are illustrated in Fig. 5.5 for two 

different pressures (P^ = 1 and 3 atm) and x-locations (x=5 and 10 cm). It 

is noted that the increase in pressure has dramatic effects on the radiative 
interaction. The conduction and radiation heat transfer results are compared 
in Fig. 6 for P = 3 atm and for two different x-locations (x * 5 and 10 cm). 
The results demonstrate that the conduction heat transfer is restricted to the 
region near the boundaries and does not change significantly from one x- 
location to another. The radiative interaction, however, is seen to be 
important everywhere in the channel, and this can have significant influence 
on the entire flowfield. The results presented in Figs. 5. 4-5. 6 should be 
physically symmetric; but, due to the predictor-corrector procedure used in 
the McCormack's scheme, they exhibit some unsymmetrical behavior. 

For the parallel plate geometry, a comparison of the divergence of 
radiative flux for gray and nongray models is presented in Fig. 5.7 for two 
different y-locations (y = 0.2 and 1.6 cm). The physical and inflow 
conditions in this case are exactly the same as for Fig. 5.3. The gray gas 
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formulation is based on the Planck mean absorption coefficient which accounts 
for the detailed information on different molecular bands. As such, this 
approach is referred to as the "pseudo gray gas formulation." As mentioned 
before, Tien and Lowder's correlation (Ref. 6) is used in the nongray formula- 
tion. It is noted that the results for the gray model are about 10-20 percent 

higher than for the nongray model. For the physical conditions of the 

problem, no significant difference in results is observed for the two y- 

locations. Both gray and nongray results are seen to increase with x because 
the pressure and temperature, in general, increase in the flow direction. The 
solution of the gray formulation in ODE form proves to be about ten times more 
efficient than the solution of the nongray formulation on the VPS-32 computer 
(the gray formulation uses 0.056 CRU's per iteration while the nongray 
formulation uses 0.57 CRU's per iteration). As such, all other results 

presented in this study have been obtained by using the pseudo gray gas 
formulation. 


The second geometry (Fig. 5.2) was selected to study the effects of 
shocks on the radiative heat transfer. The physical dimensions considered for 
obtaining specific results are L * 3 cm, Xj ■ 3 cm, X 2 * 3 cm, L x = 10 cm, and 

a - 10°. In general, the Inlet conditions considered are the same as for the 

parallel plate geometry, i.e., P^ ■ 1 atm, T^ * 1,700 K, f h 0 “ ^* 5, 

Results for this case at steady state are given in Figs. 5.8 and 5.9. 
For = 4.5, only the resulting pressure contours are shown in Fig. 5.8 and 
the results for the radiative heat flux are illustrated in Fig. 5.9 for 

= 3 and 4.5. The pressure contours clearly show the existence of the 

shock and expansion fan. The results presented in Fig. 5.9 show a significant 
Increase in the radiative heat flux over the ramp for = 4.5. The 

oscillatory behavior observed towards the end of the channel for M * 3 is 
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due to the effect of the shock reflection from the upper boundary. At both 
Mach numbers, only a slight variation in results is noted for different y- 
locations in the channel. 

For the physical conditions of Fig. 5.2 and M^ * 4.5, the radiative 
flux results along the channel are compared in Fig. 5.10 at two different y- 

locations (y 55 0.0 and 1.5 cm) and for two different compositions of gaseous 

mixture, 25* H 2 O + 75* air and 25* H 2 0 + 25* OH + 50* air. The results show 

that the radiative interaction is relatively higher for the case of 25* H 2 0 

than for 25* H 2 0 + 25* OH. It is quite likely that for the temperature and 
pressure range over the ramp, OH is highly effective in absorbing the 
radiative energy. The results for two different y-locations are seen to be 
essentially the same. It is possible that for the physical conditions of the 
problem, the radiative transfer process is closer to the optically thin limit; 
and, in this limit, the radiative flux is independent of the y-location (Refs. 
29 and 30). 

To investigate the effects of radiative energy transfer in chemcially 
reacting flows, premixed hydrogen and air with an equivalence ratio of unity 
was selected. Specific results were obtained again for the geometry of Fig. 
5.2 with the inlet conditions being exactly the same as for Fig. 5.8. For the 
physical conditions of the problem, the radiation participating species 
produced due to chemical reaction essentially are OH and H 2 0. The radiative 
interaction is started at about X/L x = 0.20 to make sure there are significant 
amounts of OH and H 2 0 produced by the reaction for active participation. The 
pressure contours for this case are shown in Fig. 5.11. A comparison of the 
pressure contours for the nonreacting and reacting flows as given in Figs. 5.8 
and 5.11 shows that the shock angle has increased in the case of the reacting 
flow. This is due to a relatively thicker boundary layer and changes in the 
thermophysical properties of the mixture. 
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The distribution of various species at the lower boundary and center of 
the channel are illustrated respectively in Figs 5.12 and 5.13. Due to the 
high temperature, OH increases very rapidly over a short distance from the 
inlet and then remains constant for both y-locations. The concentration of OH 
is found to be higher near the boundaries due to relatively higher tempera- 
ture. The rate of radiative flux along the channel is compared for this case 
in Fig. 5.14 for the reacting and nonreacting flows. For the nonreacting 
flow, the radiative flux is evaluated for 25* H 2 0; while for the reacting 
flow, it is evaluated for H 2 0 and OH produced by the actual reaction. It is 
noted that although the amount of participating species Is less in the 
reacting case, the radiative flux is considerably higher in this case than the 
nonreacting flow. This is due to the increase in temperature caused by the 
chemical reaction. 

The variations in temperature and species concentration along the channel 
are illustrated respectively in Figs. 5.15 and 5.16 for chemically reacting, 
and chemically reacting and radiating flows. The results are obtained for 
exactly the same conditions as used in Figs. 5.11-5.14. The temperature 
variations in Fig. 5.15 are given for different y-locations. It is noted that 
the temperature within the boundary layer (y = 0.025) is about five percent 
higher for the reacting and radiating interaction. For other locations (y = 
0.2 and 1.5 cm), no significant difference in the results for the two cases 
was noted. The variation in mass fraction for different species is shown in 
Fig. 5.16 within the boundary layer (y - 0.025 cm). It is found that due to 
the radiative interaction, the concentration of OH increases by about five 
percent and the concentration of H 2 0 decreases by the same amount. It should 
be noted that the radiative interaction has no significant effect on 0 2 and 
H 2 . 
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Fig. 5.16 Mass fraction vs. x for reacting and radiating flows. 
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6. CONCLUSIONS 

The two-dimensional spatially elliptic Navier-Stokes equations have been 
used to obtain solutions for chemically reacting and radiating supersonic flow 
between two parallel plates and in a channel with a ten-degree ramp at the 
lower boundary. The inlet conditions used for specific solutions correspond 
to particular flow conditions of a scramjet engine. For both physical 
geometries, diagnostic solutions were obtained to investigate the influence of 
the radiative interaction without considering any chemical reaction. 
Different amounts of H 2 0 and OH were used with air for parametric studies. It 
is noted that the radiative interaction increases with increasing pressure, 
temperature, species concentration, and Mach number. In the case of flow 
without chemical reaction, most of the energy transferred is by convection in 
the direction of the flow. As a result, the radiative interaction does not 
affect the flow field significantly. Some important results were obtained by 
considering the reacting flow in the channel with a ramp. The results reveal 
that, in the case of flow with strong shocks, radiation can have significant 
infuence on the entire flow field; however, the influence is stronger in the 
boundary layers. It is found that the numerical scheme based on the pseudo 
gray gas formulation for the radiative flux is highly efficient as compared to 
the scheme based on the nongray gas formulation, especially for the vector 
processing computers such as VPS-32. For further study, it is suggested that 
the present formulation be extended to include two-dimensional radiative 
interactions and provide parametric investigations for higher inlet Mach 
numbers. 
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APPENDIX A 


DERIVATION OF CONDUCTION HEAT FLUX TERMS 

To simplify the Eqs. (2.3) to Eqs. (2.4), the Lewis number is assumed to 
be unity. This simplification is carried out in detail for Eq. (2.3b) and the 
same is applied to Eq. (2.3a). Using the expressions for the thermal 
diffusivlty (a) and Lewis number (L e ), Eq. (2.3b) can be expressed as 

a « K/pC p 

L e = a/D 

— ftT ® &f $ 

q cy * " pD (L e C p ”5y + ^y“ V (A - 1} 

Defining the binary diffusion coefficient D in terms of the Pandtl and Lewis 
number Eq. (A.l) can be expresses as 


Pr = v/a 


a _ v/Pr p 


m bf. 


ii - ui j 

q cy * tC p 3y + ^ W h i 3 


where 


m 

C = Z f . C 
p 1=1 > p, 


The static enthalpy of the mixture is given by the relation 

h = E [h? + / T C n dn] f. 
i=l 1 o p i 1 


(A. 2) 


(A. 3) 


It should be noted that q is a dummy variable employed to evaluate the 
sensible enthalpy. Using the Leibnitz formula Eq. (A. 3) is differentiated to 
obtain 
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m 

= E 
1=1 


+ f , 


n T 

[(h? + / C ( ti) drj) 

1 0 P 1 

T 

! 0 d ” * f i C Pl 


af. ah? 

-5— L + f . — 1 

ay i ay 
(T) f] 


(A. 4) 


The coefficients of the first differential on the right side is equal to h i 
and the second and third terms are identical to zero, therefore, Eq. (A. 4) 
reduces to 


or 


ah 

ay 


m 

E 

1=1 


[h 


a f i 

i ST + f i 


'1 


aT 

■5y 


] 



m 

E 

i=l 




aT 

ay 


Substituting Eq. {A. 5} into Eq. (A. 2), q C y is expressed as 

n - . ^ ah 

q cy "7? ay 


or 



(A. 5) 
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APPENDIX B 


EXPRESSIONS FOR RADIATIVE FLUX 


The general equation of radiative flux is obtained from Ref. 5 as 


q R * - 2% j I x (0,p)e pdp -2 it f I x (t qX , -p)e p-dp. 


* 2 ^ e bX (t) G X (t) ] E 2 (t x ' t) dt 


T x \ A. 


(B.l) 


where the E n (t) are the exponential integral functions defined by 


E n ( t) = J p n 2 e dp 


-t/p 


For a diffuse surface the quantities I* (0,p) and I x ( T oX »i A ) are given as 


, + "V* 1 

2n f I x (0,p)e * pdp = 2B^ E 3 (x x ) 


1 

f 

o 


(B.2) 


1 . ~ 

X '*0X‘ • 


2* / ir lx,. -I*) e oX x 


^ ■ 28 2X E 3 U oX ‘ V 


Substituting Eq. (B.2) into Eq. (B.l), and neglecting the scattering terms 
(r x =0, p x = < x + y x ), there is obtained 

q RX (x X } = 2B 1X E 3 (x X> " 2B 2X E 3 * x oX ” T X ) 


+ 2 f X e bx (t) E 2 (i: x -t) dt - 2j C ° X e bx (t) E 2 (t-T x ) dt (B.3) 
0 T X 


59 


From the Appendix B of Ref. 5, one finds 

I E n (t) dt * - E n + l (t >: E „<°»=n^I 

Now, consider the first integral in Eq. (B.3) 

*1 = ^ E 2 ( V t} dt 

By defining y=T^-t, dy = - dt, Eq. (B.5) is expressed as 

Ii * / E 2 (y) ( ' dy) * / ^ E 2 (y) dy 

T X 0 

- -E 3 (y)|V E 3 (°) -E 3 ( V 

Thus 

E 3 ( V - e 3 (0) - h - t! dt 

The second integral in (B.3) is written as 

Ij> - / T ° X Ep ( t“T. ) dt 

By defining y = t-x^, Eq. (B.7) Is written as 

I, - / T ° X * X E,(y) dy = -E 3 (y)| T ° x 
0 1 0 

■ - t E 3 ( V-V - E 3 (0 »] 

Thus, 

E 3<*0 - V * E 3 (0 > -h-\- E 2«»-V 

T \ 


(B.4) 


(B.5) 


(B.6) 


(B.7) 


dt (B.8) 


A substitution of Eqs. (B.6) and (B.8) into Eq. (B.3) results in 
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< W T \ ) = e lX " e 2\ + 2 ^ e b\ (t) “ e l\^ E 2 (t \ “ tJ dt 


- I ° X Ce 5x (t) - e 2x ] E 2 (t - t x ) dt} (B.9) 


Equation (B.9) is exactly like Eq. (3.9) if it is written in terms of the wave 
number (w) instead of the wave length (X). 
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APPENDIX C 

BOUNDARY CONDITIONS FOR GRAY MODEL 

The evaluation of the second order differential equation, Eq. (3.18), 
requires two boundary conditions. In this Appendix, the boundary conditions 
for nonblack diffuse surfaces are derived. 

Applying the exponential kernal approximations to Eq. (3.8), one obtains 


d RX 










( V t) 


dt 


3 

7 



e bX (t) 


3 

7 


(t-T x ) 


dt 


(C.l) 


By employing the Leibnitz formula, the divergence of the radiative flux is 
expressed as 


dd RX ! V 
cTT 


3 R / 7 - 3 B e" 

7 B 1X e 7 B 2X e 


7 ^ox'V 


9/4/ o S b x(t) e 7 ♦iw 


Q x - * 4 (t-T x ) ? 

^ / e b x(t) © dt + <ir e K ^ (*c % ) 

T, 


7 bx v V 


(C.2) 


Evaluating the radiative flux terms at the lower and upper walls the relations 
for q RX and divq RX are expressed as 

- 3 x x - 3 t 

«RX (( -' * B 1X * B 2X e ^ °" - 1 V (t > e ^ dt (C - 3) 


* 
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< W T oX ) = B 1X e ‘ B 2X + 7 J ° X e bX (t) e 


7 ( ' c oX" t) 


dt (C.4) 




3 a 3„ " 7 T oX 9 /ox a 

7 B 1X 7 B 2X e ‘ ? J e bX 

0 


dt + 3 e b!l (0) 


(C.5) 


dq R ( V, . 3 „ ’ 

"J%1, r •'7 B U e 
x W 


7 T oX 3 n 9 
" 7 B 2X “ T 


b - 9 f T oX * m “ 7 (T ox _t) 

B 9\ 7 j e 


+ 3 V'V* 


(C.6) 


Multiplying Eq. (C.3) by | and Eq. (C.4) by (- |) and then substituting 
into Eqs. (C.5) and (C.6) respectively, one obtains 


1D1 , O 

rrl + 7 q R\ {0) = 3 B u " 3 e bx (0) 


X 1 T^=0 


(C.7 ) 


x 1 W 


*7 W T oX> = 3 e b\ (T ox> * 3 B ; 


(C.8) 


The radiosities B u and B 2x are expressed as 


B 1X = e lX e blX + ^ e i/ ~ q D1 (0)] 


IX' L IX h RX v 


(C.9a) 


B 2X - e 2X e b2X + ^ - Qdi 


2X 2X H RX x W 


(C.9b) 


Rearranging Eqs. (C.9) and substituting for B u and $ 2X into Eqs. (C.7) and 
(C.8) the boundary conditions are expressed as 


( 7T7 " 7 } q RX {0) " 


1 d ^ RX | 

1 V° 


e bix - e bx (0) 


(C.lOa) 
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' 7 } q RX 


+ 


1 dq RXi 


V T oX 


e bX (T 0\ ) 


e b2X 


(C.lOb) 


For the gray medium, X dependence is deleted; and then Eqs. (C.10) are 
exactly like Eqs. (3.19). 


